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Abstract 

The derivation of the Navier-Stokes equation starting from the Liouville equation using projection 
techniques yields a friction term which is nonhnear in the velocity. Using the results of multilinear 
mode-coupling technique for calculating equlibrium correlation functions, it is shown that the second- 

P3 ' order part of the term is non-zero, thus leading to an incorrect formula for the equation. 

li I Key words: statistical thermodynamics, projection operator technique, multilinear mode-coupling 

OJj, theory, hydrodynamic equations 
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o , 
'Si 1 Introduction 

i-Q , The derivation of hydrodynamic equations by Zwanzig-Mori projection technique is a well-established 

^~j- method of statistical mechanics; see, e. g., the presentations in [T], [2]. But analyzing the results further 

leads to a problem which is connected to the fact that the Navier-Stokes equation is of second order in the 

velocity. By derivation, the friction term is nonlinear in the velocity; in order to keep the equation correct 

^ I it would be necessary for the second-order part of the term to vanish. But, as is shown in this paper, 

the derivation furnishes an expression which is generally non-zero. Since the Navier-Stokes equation is 

?C^' ] obviously a correct formula, the question arises which features of the derivation process provide these 

PsJ ■ undesired details. 



o 
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c^ 



2 Some basic formulas 



In this section some definitions and formulas are listed which are found in the textbooks of statistical 
mechanics. We are working with a simple fluid, i. e. we are looking at a system of N point particles 
^^ . with positions y^ and velocities v^ which collect to the phase space variable z. Indices which number 

?—( I the particles are denoted by j, k,- ■ ■ . The space densities of the conserved variables mass, energy and 

momentum, called n, e, p, are collected to a matrix a: 

N 

aa(x) =^aja(z)(5(x-yj) (2.1) 

'y^' + ^E^dy^-y'i)'"^^^-! (2.2) 

The greek index a runs over 1, 2, a where a runs from 1 to 3 so that a attains 5 scalar values in total. 
m is the particle mass, if{\yj — y/|) the intermolecular potential which I write fji for short. The motion 
of the particle system is described by the variable z as a function of time t . From this, all phase space 
functions, especially the variables a, are functions of time described by the Liouville equation: 



f^£« (2.3) 



^=^T.'^^^,H-yy^-yy^H-V.^) (2.4) 

The Liouville Operator C is presented here as a real differential operator. H is the Hamilton function: 






(2.5) 



The formal solution of the Liouville equation reads: 

a(i) = e*^* a(0) (2.6) 

The conservation property of the a is expressed by: 

^ = -V • s = -Vcs„, (2.7) 

s are the flux densities. Paired indices are summed over. The s show a formal structure similar to 
the a; cf. ()2.ip . The momentum flux density reads: 

Sa6(x) = ^ I mV.aV,,, " 2 ^ ^ 'I'-^'l ^^^^^ J ^^"^ ~ ^ '^ ^^'^^ 

Here r^/ = yj — y/ . - The expectation of a phase space function A{z^ N) with respect to the probability 
density /(z, N) is defined: 

(A)=ir{A(z,iV)/(z,iV)} (2.9) 

The definition of the Operation tr{} reads generally: 

tr{A{z,N)}^'^ dzA{z,N) (2.10) 

N=l'' 

3 Hydrodynamic equations 

For the projection calculations the formalism in [1] is adopted where the time dependence of the dynamic 
variables is considered while the expectiations are defined with respect to the probability density at time 
t = (Heisenberg picture). Moreover, in exponential operators where the projector V and the Operator 
C appear together in the exponent, C appears 'before' (left of) V. In this paper, this is called the C-V 
formulation which stems from Zwanzig (see |31 where the Schrodinger picture is used; the succession of 
Operators reverses if one switches to the Heisenberg picture) while Mori 31 seems to prefer the V-C 
formulation. 

The result of the derivation in [T are the formulas (8.1.13) together with (8.1.12), (8.1.7) there; in 
the denotation here: 

^^ = -^c{s^c)L,t + D^it) (3.1) 

The expression for the dissipative force reads: 

-D„(x, t) = -Vc f dt' I dx'i?„0ed(x, i, x', i')Vd&/3(x', t') (3.2) 

with the kernel function: 

Rc,pcd{^, t, X', t') = ([g{t', t)s„e(x, t)]s^d(x', t'))L,t' (3.3) 

Oi.t denotes the expectation with respect to the probability density of the local equilibrium fhit)'- 



hit) - V(A^) exp(vl/(i) - a(z) * b{t)) (3.4) 

where 

/i is Planck's constant. * denotes an operation which consists of a product, a summation over a greek 
index and a space integration; the respective variables and indices are not written down. The conjugated 
parameters bp{yi,t) are: 

&-{/3(fw'-/^),/3,-/3u} (3.6) 

here (3 — l/{kBT); ks is Boltzmann's constant and T the temperature; these as well as the chemical 
potential /i und the fluid velocity u are generally functions of x and t. ^{t) is a normative quantity which 
ensures tr{/L(i)} = 1. In p.3p . s denotes the reduced flux densities: 

s(x,t) = (l-7'(i))s(x) (3.7) 

V is the Zwanzig-Mori projection operator; for any phase space funktion g, it is defined: 

Vg^{g)L + {9 5a)L*{5a5a)l^*5a (3.8) 

5a — a — (a)L- ()^^ denotes the inverse matrix. Generally ()i, and therefore V , depend on time. 
We have R — for a — I or/and [3 = 1, since then the reduced fluxes are zero. G{t' , t) is a time-ordered 
exponential operator: 

g{t', t) = exp_{ / dt" C{\ - V{t"))} (3.9) 

We need formula ([S3]) for a = l,a. One obtains ([Tl (8.4.1), (8.3.12) und (8.4.12)): 

{sic)L,t^ puc (3.10) 

{Sac)L,t ^ PUaUc + PSac (3-11) 

p is the macroscopic mass density, and P the pressure. - From here on the processes considered 
are restricted. The Problem mentioned in the introduction appears even under this specialization, and 
the necessary formulas are considerably simplified. Constant mass density/temperature processes are 
considered; then, the chemical potential /i ist constant too. From (13. 6|) . the term /3 = 2 in ()3.2p is zero; 
thus, the sum runs over /3 = 6 only. Finally, I confine to stationary currents, i. e. u =:const(i). Then, 
expectations of single phase space functions as well as the conjugated parameters and V are constant in 
time; G reduces to a non-ordered exponential operator. Finally, the upper limit of the time integral in 
p.2p may be extended to infinity. Then it is reasonable to take the time integral under the definition of 
the kernel function. For a = a, one obtains, instead of p.2p . p.3p : 

Da{x, t)=\/c f d'x Rabcdi^, x')V>b(x') (3.12) 

Rabcdi^, x') = /3 / dt(e^(i-^)* sac(x)] Jm(x'))l (3.13) 

Finally, from (|3.ip . together with p.lOp . (|3.1ip . one obtains for stationary constant density/temperature 
currents 

V • u = (3.14) 

pUc\/cUa = -^aP + Da (3.15) 

with p.l2p . (|3.13p for D. These are the incompressibility condition and the momentum equation for 
stationary constant density and temperature current. 

R and D are nonlinear functionals of the velocity field. From (j3.12p it is obvious that in lowest order 
D is linear in u. One obtains this order when one takes i? at u = 0. u enters R exclusively via the 
formula for the local equilibrium (j3.4p . We have: 

fLiz,N)\^^o = foiz,N) (3.16) 



/o is the probability density of the thermodynamic equihbrium for prescribed mass density and tem- 
perature: 

/o(z, N) = ij{N) exp(*o + |3{^^N - H{z))) (3.17) 

^0 is the norming quantity corresponding to ^. i? at u = is called R^^^\ 

4°L(x-x') =/? / di(e^(i-^«)*(^o)ac(x)](^o)b.(x')>o (3.18) 

Jo 

(So)ac = (1 - Po)Sac (3.19) 



'Poa = {9)0 + {9 5oa)o * (SoaSoa)^ ^ * Squ (3.20) 

()o is the expectation with respect to total equilibrium, and we have Soa = a— (a)o. Since correlations 
in total equilibrium are transitive, i?^"^ depends on the single space variable x — x' only. Provided the 
space integral over i?'^"^ exists, formula (|3.12p . in the linear approximation, can be localized. The space 
integral resembles the quantity in [T], Gl. (8.5.20); though the formula is obtained there in a somewhat 
different way. The integral obeys several symmetries, which lead to the expression jT] (8.5.21); in the 
denotation used here (D^^^ is the linear part of D): 

^i'^(x)=7abcdVcVdUb(x) (3.21) 

labcd = / d^Rabidi'^) = i^abScd + 5 adhc^ + 5 achd{^ - g??) (3.22) 

77, c are the (dynamical) shear and bulk viscosity. If this is introduced into (I3.2ip . and (|3.14p is allowed 
for, 1^ eliminates, and ()3.15p changes to: 

pUcVclia == -VaP + ll^'^Ua (3.23) 

This is the Navier-Stokes equation fiir stationary constant density/temperature flow. 

4 2^'^ order term of the friction force 

The Navier-Stokes equation p.23p is of second order in the velocity. In order that it results correctly 
from the momentum equation ()3.15|) as an approximation for small Reynolds numbers, the 2"'' order part 
D'^'^^ of the friction force must vanish. This part is built in p.l2p with the linear part R'^^^ of the kernel 
function (I3.13p . By Taylor's theorem, the latter reads: 



i?ilL(x,x') . |dx"^^^|;;i^|„=0Ue(x") 



(4.1) 



It is useful to express the functional derivative with respect to u by the derivatives with respect to 
the conjugated parameters 6^ p.6p . These are parametrically related to ; since 6i is of 2"'' order in u 
and 62 — const(u), for u = only the index value e = e remains, and we have: 

^i?a6cd(x,x') _ 5i?Qbcd(x,x')| , . „. 

5uM") '"=° "^ <5&e(x") '"=° ^^-^^ 

The calculation of the functional derivative is performed in the appendix. The result is (|A.15P : 

(5i?afccd(x,x') 



<56e(x") '"- 
With gU, g3), this yields: 



P / dt([e(i-^°)^*(so)ac(x)](Jo)6d(x')Pe(x"))o (4.3) 



4L(X,X') = -/3/rf-"^^||^lu^0Ue(x") 

= /J^y dx" ('y"dt([e(l-^«)^*(so)ac(x)](so)M(x')Pe(x"))o) ^e(x") (4.4) 



Therefore, D^^^^ reads: 

= /jV.y dx'l dx" ('y"di([e(i"^'')^*(so)„c(x)](Jo)bd(x')Pe(x")>o) u,i^")V,Uh{^') (4.5) 

4.1 Calculation of the kernel function 

In this subsection, the calculation of the kernel function in (|4.5p (the quantity in parentheses) is performed, 
using a result of multilinear mode-coupling theory. Instead of (|2.ip . orthonormal phase space densities 
are used. By replacing the energy density e by a suitable linear combination of e and n, the densities can 
be made orthogonal: 



d-x{Soaaix)Soa*p{0))o = Xa^ais (4.6) 

The orthonormal variables ha read: 

ha = —= (4.7) 

\ Xa 

The flux densities conjugated to the h^ are called Vac, the corresponding reducted quantities fac = 
(1 — Vo)rac- By introducing the new definitions into (j4.5p one obtains: 

= /3^piVey dx'l dx" ('y°°dt([e(i-^°)"f,,(x)]fM(x')Pe(x"))o) ^e(x")V>fc(x') (4.8) 
The calculation is done in Fourer space. We have: 

i^i^) (k) = /3 ^f _i_ J dqNabcO^, q, k - q)ubiqW{k - q) (4.9) 

iVafcc(k,q,q')-fcdge / dt([e(l-^o)£*^^^(k)]^*^(q)/i*(q'))o (4.10) 



The kernel function N is to be calculated. The multilinear mode-coupling technique is for calculating 
correlation functions which are defined with the original exponential operator expjLt}, that is, without 
1 — Vq. In order to find a connection, one uses the operator identity: 

Jo 
Application results in the formula: 

Nabc{k,q,c^)^iqeJ dtlj dt'f„,(k,t-i')([e^*'^.(k)]rL(q)/i:(q')) -«fcd([e''*'w(k)]r,*Jq);z:(q')) 

(4.12) 

f„^(k,t) = Kkd{^ac{tWpd){^)W=0 (4.13) 

l>„,(k,i)=e(i-P«)^*f„,(k) (4.14) 

{AB*){k) is defined to be the Fourier transform of the correlation (yl(x)_B(0)). Tap is the localized 
memory function of the process. The ~ points to the fact that the function $ is defined with the aid of 
the projection operator Vq which is linear in the microscopic densities ha (see below for the discussion of 
the projection operator of the multilinear technique). It is assumed that the time integral of the memory 
function exists: 

/■CX3 

7„^(k)=/ dtfafAKt) (4.15) 

Jo 



Then, f in (|4.12p can be "localized" in time. 7^^ is found to be diagonal; thus, ii a — a, then /3 — b. 
The other correlations appearing in (|4.12p can be related to the triple correlation C3 of the orthonormal 
densities: 

{C3)c,p^{t) = {[e"^' h^]h}ih;) (4.16) 

Finally, for the application in (|4.9p it is possible to replace Nabc by the symmetrized form ■^{Natc + 
Nacb) which is denoted by the same symbol. One obtains from (|4.12p : 



"^Nabc = -fiae {C3)ebci0)- dt{i UJbT{C3,)eTc{t) + i l^ cT{C3,)ebT{t)) 

Jo 

-^^^^^^\t=0-^iObr{C3)arM-^iU,riC3)abr{0) (4.17) 

at 

Kaf3 ^Ia0+'>'^al3 (4-18) 

w„/3(k) = kdivadh}) {k)\k=a (4.19) 

N and C3 depend on three wave numbers fc, q, q' . The wave numbers of the other quantities in 
(|4.17P are7^,(k), Wae(k), Wf,^(q), Wcr(q')- 

By (|4.17p . we connected the kernel function N to the triple correlation C3 which is defined with the 
exponential Operator exp{Lt}. The formula for the latter we take from the multilinear mode coupling 
theory for correlation functions |B]. This theory is worked out in Fourier space before performing the 
thermodynamic limit. Then, the correlations are proportional to the Volume V. In what follows all 
correlations are divided by V; then, V does not appear any more in the formulas. From the formulas 
(31), (30) there, one obtains, with the denotation used here: 

(^3)0/37(0 = iC2)Ut)JefJy - I dt'{C2)Ut - t')lkdS,prd{C2)p0{t'){C2)r-,{t') (4.20) 

JO 

{C2)c.f}{t) ^ {[e"' K]h*p) (4.21) 

Sap-fd = (»'Qrf^^/l^)|q,k,k'=0 (4.22) 

Ja(3-i — (C3)Q,/37(0)|q^k,k'=0 (4.23) 

For C2, one obtains, after "localization" in time corresponding to (|4.15p : 

(^2)0/3(0 = -Hc.^{C2hp{t) (4.24) 

Ka/3 = 7q/3 + « ^a/3 (4-25) 

7 is the dissipation matrix defined corresponding to ()4.15p but with the projection operator of the 
multilinear Theory which is nonlinear in the ha- When one introduces (I4.20p into ()4.17p , one obtains 
an expression for N which contains the matrices w, 7, 7, J and S . For the present intention to check 
whether the 2"'' order part of the friction force is different from zero, one can neglect the difference 
between 7 and 7. Then one finds: 

Nabc = l^ilbp^ca + SbpJccr)il^peSar + 5 peK^rY^ ikdSaerd (4.26) 



The kernel function N will be used to calculate the friction force term D2 (|4.9p . with the incompress- 
ibility condition p.l4p k • u = to be incorporated. This being accounted for, for N there remains a 
relevant part: 

Nabc = -^ikdSabcd (4.27) 



One finds: 



Sabcd = —;=={5ab6cd + SacSbd — ^Sad^bc) (4.28) 

A=^|P(^|P)- (4.29) 



and by introduction into 

I?i'^ (k) = j£r^ J rfq (ikd Ua{q)ud{k - q) - ^ika Ufc(q)wb(k - q)") (4.30) 

It is seen that the 2"'' order part of the friction force is calculated to be not zero. The formula will 
be evaluated in the next section. 

5 Consequences 

The consequences of the result (|4.30p will show up clearly if one uses the solenoidal form of the Navier- 
Stokes equation, that is, the equation with the pressure term elimited. Instead of (|3.23p . the instationary 
form of the equation is chosen, and the 2"*^ order term of the friction force D'^^\ — VcRac ist added. 
Here Rac with two indices is the stress tensor which should not be confounded with the kernel function 
Rabcd- The calculation is performed in Fourier space. The Navier-Stokcs equation supplemented by D^^> 
reads: 

p(^ + J^k^ua) = -ikaP - tk, fj£^ I dqua{<l)u,{k - q) - R^A (5.1) 

After a short calculation (see, e. g., [5] appendix D, (D34)) the solenoidal form of the equation is 
obtained: 

p{^ + ^k^ua) = -ikcSab (j£y^ J dciut{<l)Uc{k - q) - R^A (5.2) 

£ab = Sab ~ kakb (5.3) 

k is the unit vector attached to k. Now, for D2 the result (|4.30p is introduced. Since we have 
kc^ac — 0, the 2'"^ term does not contribute. One obtains: 



ikcSab i 7^-^ / (iub(q)wc(k - q) - i?^ 



(5.4) 



The 2"^^ order part of the friction force cancels the convolution term of the Navier-Stokes equation, 
so that the equation reduces to its linar part. Since by phenomenological evidence the convolution term 
is an essential part of the equation, this theoretical result cannot be correct. 

6 Summary 

The derivation of hydrodynamic equations by Zwanzig-Mori projection technique has been reviewed. 
For simplicity, incompressible constant density/temperature fluids have been considered. The friction 
force term in the momentum equation is a nonlinear functional of the fluid velocity; when the linear 
approximation is taken, the Navier-Stokes equation is obtained. In the present paper, the second-order 
term of the friction force has been calculated. It contains a three-point time correlation function which 
has been evaluated using a result of multilinear mode-coupling theory; the final result is ()4.30p . Since the 
Navier-Stokes equation is second order in the velocity, in order to obtain it properly as an approximation 
of the momentum equation for small Reynolds numbers, one would expect the second-order part of the 
friction force to vanish. Thus, (|4.30p is clearly an undesired result; it would be very valuable to know 
which detail of the derivation may be responsible for it. 

A Appendix: Calculation of the functional derivative 

In order to calculate the derivative in (|4.2p . the formula for R p.l3p is written in detail: 

/>oo 

i?abcd(x,x') =/3 / dttr{/i[e^(i-^)*(l-7')sac(x)](l-7')sbrf(x')} 
Jo 



(A.l) 



The expression depends on be 4-fold, namely, in the formula for the local equilibrium /i, and in V 
appearing 3-fold. For abbreviation, the formula is written: 






SR 
lb 



(A.2) 



In the consecutive formulas, in the first row the definition of the several parts is expressed. For the 
calculation some auxiliary theorems for projection operators are used, which can be found in the text 
books. For the 1"* partial term: 



The 2"'' term 



SRa^ 



bed 



<56e(x") 



(1) 



/•OO 

= -(3 dt([c^(i-^)*s„,(x)]sbd(x')<5pe(x"))L 
Jo 



(A.3) 



SRa 



bed 



1(2) 



<56e(x") 



= -/3^°'dt([c^(i-^)*5„,(x)]-^^^.,,(x'))L 

= f3 dt([e^(i-^)*s,,(x)]7's,d(x')<5pe(x")>L 
Jo 



(A.4) 



The Operator V acts on everything on the right (provided the action is not limited by parentheses). 
The 3'''' term: 



SRa^ 



bed 



Sb,{^" 



(3) 



-/3^°°di([c^(l-^)*^^^.„e(x)]^,,(x'))L 

/•OO 

(3 / dt([e^(i-^)*7's„e(x)<5pe(x")]sM(x'))L 
Jo 

P / di([e(l-^)^*(l - r)rSaci^)6pe{^")]sMi^'))L = 



The 4*'' term is defined: 



SRa 



bed 



(5fce(x") 



(4) .OO g^C{l~V)t 



The formula for the derivative of the exponential operator reads: 

s:pC(l-V)t ft 

^ = di'e^d-^)* /:7'<5pe(x")(l-7')e^(i-^)(*-*) 



Sbei^") Jo 



We use the identity: 



dt dt'A{t',t-t')= dt dt'A{t',t) 
Jo Jo Jo 



(A.5) 



(A.6) 



(A.7) 



(A., 



One obtains: 



SRa 



bed 



<56e(x") 



(4) 



OO /'OO 



= /?/ dt dt'([e^(i-^)*'£Wpe(x")(l-P)e^(i-^)*s„e(x)]s,,(x'))L 
^0 Jo 



OO /-OO 



^0 

/•OO /"OO 

= P dt dt 

Jo Jo 

The integration over t' is performed: 



(3 I dt I dt'([e(i-^)^*'(l-7')/:Wpe(x")e(i-^)^*s„,(x)]sbd(x'))L 

'([A e(i-^)^*' VSpei^") e(i-^)^* s,,(x)]5,,(x'))l (A.9) 



r5i?„ 



bed 



6K{^") 



(4) 



/3 / dt I hm ([e(i-^)^*' T'-^pelx") e^-^)^* s,,(x)]s,d(x'))L 

Jo ^t'^oo 

- ([7'(5pe(x")e(i-^)^*s„e(x)]sb,(x')>L} (A.IO) 



In the 1^^ term, the operation V is performed by 



( / dt lim ([e(i-^)^*>,5p,(x")e 



//^„(l-•p)£t 



Sac(x)]s6d(x'))i 



= /? 



dt hm < 

t' — >-oo 



+ {Sp,i^")[e^^-''^^'Saci^)]5a)L{SaSa)l^Sa}] sm(x'))l 
= Z?!" di {(<5pe(x") e(i-^)^* s,e(x)>i(Jbrf(x')>L 

+ ((5pe(x")[e'^-^'^*Sac(x)]5a)L(Ja(5a)^i lim ([e^"^)^*' Hs6d(x')>L| (A.ll) 

t' — 'OO J 

Here, the 2"'* factor of the 1^* term vanishes, therefore the f** term vanishes in total. All expressions 
are finely to be taken at 6 = 6o- Therefore, the limit limt'_>oo of a correlation function in thermodynamic 
equilibrium must be performed. One assumes that this equates to the product of the expectations of the 
factor functions. This assumption has been checked with the aid of the multilinear mode coupling theory, 
and verified for the special cases appearing here: 



lim ([e(i-^«)^*' Soa]{l ~ Vo)sbd{^'))o = {Soa)o{{l - ^o)sm(x')>o - 



(A.12) 



with Va by p.20p . Therefore, the 2"*^ term in (|A.11[) is zero also, so this is true for the f** term in 
(|A.10[) in total. Finally, the 2"'^ term in (jA.lOp vanishes, since one factor starts with V, the other with 
1 — v. One obtains: 



SRabcd 



(4) 



Thus, we finally find that the derivative (|A.2p . after letting u = 0, will be equal to the sum of the 
first two terms of the r. h. s.; with (|A.3p . (|A.4p : 



SRa 



bed 



5K{^") 



(1) 



5Ra. 



bed 



<56e(x") 



(2) 



(3 / dt([e^(i"^'* s„,(x)](l - V)hd{^)5pe{^'))i 



-13 / dt([e(i-^)^*s,,(x)],Sbd(x')<5pe(x"))L 
Jo 



(A.14) 



One sees that the shortest expression results with the T'-before-ZI formulation. Finally, by taking 
u = 0: 

5i?,,,,(x,x'), ^^ ^ _^ /°°dt([e(i-^»)^*(Jo).c(x)](Jo)M(x')Pe(x")>o 



5be{yi") Jo 

The last factor pe does not show a S any more since (pe)o = 0. - The calculation of the derivative 
has been performed in T'-before-ZI formulation also. The in-between steps look somewhat different; for 
instance, we have [Y^^ ^ 0. The final result coincedes with (|A.14|) . 
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